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Abstract. Axions with mass m a > 0.7 eV are excluded by cosmological precision data be- 
cause they provide too much hot dark matter. While for m a > 20 eV the a — > 27 lifetime 
drops below the age of the universe, we show that the cosmological exclusion range can 
be extended to 0.7 eV < m a < 300 keV, primarily by the cosmic deuterium abundance: ax- 
ion decays would strongly modify the baryon-to-photon ratio at BBN relative to the one 
at CMB decoupling. Additional arguments include neutrino dilution relative to photons by 
axion decays and spectral CMB distortions. Our new cosmological constraints complement 
stellar-evolution and laboratory bounds. 
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1 Introduction 

The era of precision cosmology has opened new opportunities to use the universe as a particle- 
physics laboratory. Cosmological neutrino mass limits and the prospect of eventually mea- 
suring this fundamental parameter in the sky is, of course, the most notable example [1, 2]. 
Analogous arguments allow to derive limits on other hypothetical low-mass particles and 
here axions are perhaps the most interesting case [3-7], where cosmological data imply 
m a < 0.7 eV [8]. Of course, if axions are the cold dark matter of the universe [9], they 
are very light {m a < 10 fieV) and such arguments are moot. Therefore, the hot dark matter 
limits are primarily useful to complement stellar energy-loss arguments [10] and the direct 
search for solar axions by CAST at CERN [11] and the Tokyo axion helioscope [12]. 

Axions decay by a — > 2j, the lifetime dropping below the age of the universe for 
m a > 20 eV. For a sufficiently large mass, the decay happens early enough to thermalize 
the secondary photons and hot-dark matter or overclosure arguments no longer apply. 

However, decaying axions are still constrained by several cosmological arguments. Be- 
fore the cosmological precision era, such questions were discussed in great detail for axion-like 
particles [13, 14]. A central argument was the impact of excess radiation provided by axions 
during big-bang nucleosynthesis (BBN). Today it has become clear that BBN does not credi- 
bly exclude excess radiation provided by one additional thermal degree of freedom and some 
amount of excess radiation may even be favoured [15, 16]. Moreover, cosmological precision 
data somewhat prefer excess radiation equivalent to 1-2 neutrino species at the epoch of 
CMB decoupling [17-21]. Therefore, we here address the question which axion mass range 
beyond 0.7 eV is actually excluded by cosmological arguments alone. 
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If axions decay not too early, spectral distortions of the CMB provide restrictive con- 
straints [14]. For earlier decays (larger masses), the decay photons thermalize and heat the 
CMB relative to cosmic neutrinos. Since all cosmological quantities are normalized to the 
measured CMB temperature, this means that effectively neutrinos are diluted. If axions de- 
cay non-relativistically and thus out of equilibrium, this effect can be large and leads to a 
lowered radiation density at CMB decoupling, increasing the tension with precision data that 
actually favour a radiation excess. Most importantly, the baryon density at CMB decoupling 
is also diluted. Therefore, the baryon abundance relevant for BBN is larger than implied 
by CMB data, leading to a smaller primordial deuterium abundance. Confronting this with 
observations requires that axions are not significantly abundant during BBN. This new BBN 
limit requires m a > 300 keV. Our new bounds therefore nicely overlap with the constraints 
from beam dump and reactor experiments [22-27]. 

To derive this result we first describe the relevant phenomenological properties of axions 
in Sec. 2. In Sec. 3, we discuss the modification of the cosmic expansion history and the rela- 
tive contributions of different forms of radiation in our decaying axion scenario. Constraints 
based on BBN, neutrino dilution, and CMB distortions are presented in Sees. 4-6 before 
concluding in Sec. 7. 



2 Axion properties 

The axion arises in the Peccei-Quinn solution of the strong CP problem as the Nambu- 
Goldstone boson of a new global spontaneously broken U(1)pq axial symmetry [28]. The 
main ingredient of the PQ mechanism is that U(1)pq is colour anomalous, providing for an 
axion-gluon interaction of the form 

C-aGG = TT" G% u G a n V . (2.1) 

J a OTT 

The energy scale f a , the axion decay constant, is the main parameter that fixes low-energy 
axion phenomenology. This term automatically cancels the CP-violating 6 term of QCD 
because it generates nonperturbatively a potential for the axion, driving the axion field to 
the CP-conserving position. The axion mass thus created is 

m „ = &p^ 6eV (~), (2.2) 

fa l + Z V fa ) 

where m n = 135 MeV is the neutral pion mass, f n = 92 MeV its decay constant, and 
z = m u /md. We use the canonical value z = 0.56, although the possible range is z = 0.35- 
0.6 [29], implying a 10% uncertainty of the axion mass. 

In the minimal scenario we shall consider, axions do not have tree-level interactions with 
leptons. There are thus called hadronic axions. The freeze-out of sub-eV hadronic axions is 
determined by their interaction with pions. However, we here consider much larger masses 
and much later freeze-out. Their only relevant interaction is with photons, 

£ a77 = i^F^F^a = -g aj E B a , (2.3) 

where F is the electromagnetic field-strength tensor, F its dual, and E and B the electric 
and magnetic field, respectively. The coupling constant 

a ( E 2A + z\ a 1 n . (n A . 

^ = ^tAn--,—,) s ^ s ' (2 ' 4) 
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consists of a model-independent contribution from a-ir-rj mixing and a model-dependent one 
parameterized by the ratio E/N of the U(1)pq electromagnetic and colour anomalies. The 
uncertainties from the up/down quark- mass ratio z and from the model-dependence of E/N 
can be lumped into the parameter 5. We have normalized 5 such that 5 = 1 for E = 0, 
happening for instance in the KSVZ model [30, 31]. 
The well-known a — > 2j decay rate is 

r — = 4^ * Llx10 "" ■"'(!£) v - (2 - 5) 

For 5 = 1, axions with m a < 18 eV live longer than the age of the universe. In this paper, 
we shall consider 5 values of order one, excluding radical fine-tuning or very large E/N 
scenarios. Note however, than a value 5 = provides the most restrictive limit: axions would 
not disappear and eventually over-close the universe. 

In extended models such as the DFSZ [32, 33], axions also have tree-level interactions 
with leptons. The Nambu-Goldstone boson nature of axions determines the form of the 
interaction to be of the derivative axial-current form [34] 

£l = irf ^Z7 M 75^Z (dpd) = ^-p- tpnstpl a, (2.6) 

^Ja J a 

where the Cj's are model-dependent numerical coefficients . Note that the coupling is propor- 
tional to the lepton mass 1 . For late cosmology we will be mostly interested in the coupling 
to electrons since the coupling to neutrinos is extremely suppressed by the small neutrino 
masses. 

3 Decaying axion cosmology 

Our main cosmological arguments depend on modified baryon-to-photon and neutrino-to- 
photon ratios at various epochs. One difficulty is that the freeze-out epoch of axions in the 
relevant mass range can coincide with the e + e~ annihilation epoch and with the BBN epoch 
itself. Moreover, for m a < 20 keV, axions first decouple from the Primakoff process and then 
recouple by inverse decays, causing a somewhat complicated evolution of the radiation fields. 
As a first step we therefore discuss this evolution. 

Thermal axions are produced after the QCD epoch by hadronic processes such as tt+tt f-> 
7r + a. As the universe cools, hadrons disappear and the photo-production of axions from 
electrons, e ± 7 —> e^a, dominate. For hadronic axions, the most relevant reaction is the 
Primakoff process (see Fig. 8a in Appendix A). A simple estimate for T S> m a is 

T P ~ ag 2 ai n ei (3.1) 

where n e is the density of electrons plus positrons. For non purely hadronic axions, the direct 
axion-electron coupling allows for the Compton process (see Fig. 8b in Appendix A). An 
estimate of its rate is 2 

T C ~ a n e . (3.2) 



1 The equivalence between the axial and pseudoscalar forms of Eq. (2.6) is not always guaranteed but it 
will be so for the applications of this paper where we consider processes involving the emission or absorption 
of one Nambu-Goldstone boson at a time. 
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At temperatures close to the electron mass, this process is more efficient than the Primakoff 
since involves less powers of a, and when present will dominate the axion thermalization. 
Other processes like e + e~ — > ja are subdominant and can be neglected in the following. 
More details about the Primakoff and Compton processes are given in Appendix A. 

The Primakoff and/or Compton processes freezes out when their rate drops below the 
cosmic expansion rate 

u _ldR_ I 87T _ I 8tt tt 2 2 

where R is the cosmic scale factor, mp\ = 1.22 x 10 19 GeV is the Planck mass, p the total 
energy density of the universe, T the photon temperature and the last equality defines g*, 
the effective number of thermally excited degrees of freedom. The rates Tp and Tq decrease 
exponentially when e + e~ annihilation sets in, so the Primakoff and Compton processes indeed 
freeze out (see Fig. 1). 

Another processes based on the two-photon vertex are decay and inverse decay. For 
T S> m a the inverse decay rate is 

m„ — Ami i m n \ 

77— >a — - 1 - a— >77 o \ / ; V° ,4: / 

where {m a /uj) is a thermally averaged time dilatation factor. Decay or inverse decay is 
only possible if m a > 2m 7 , where m 7 is the plasmon mass. For T S> m e , we have m 7 = 
\J?>a'K j2 T < T while for T <C m e we have m 7 <C T so the decay/inverse decay channels 
open up not far from m a ~ max{T, m e }. Leaving aside the threshold, for relativistic axions 
r 77 _> a oc T _1 and this process actually recouples at a certain temperature (thin solid line in 
Fig. 1). 

In Fig. 1 we show the decoupling and recoupling temperatures (horizontal axis) as a 
function of axion mass (vertical axis). For a given m a the cooling universe moves on a hori- 
zontal line from right to left, beginning with axions in thermal equilibrium (yellow shaded). 
For m a > 20 keV the horizontal line never leaves yellow territory and axions remain in ther- 
mal equilibrium forever, first by the Primakoff process and later by inverse decays. If axions 
couple to electrons, this region is extended until m a ~ 10 keV through the Compton process. 
In this case, the cosmic axion history is governed by thermal equilibrium. The populations 
of all forms of radiation are determined by entropy conservation (adiabatic regime). 

Another relatively simple range is m a < 1 keV where Primakoff freezes out when the 
e + e~ plasma is fully populated and Compton is not very effective. The approximate temper- 
ature range of e + e~ annihilation is delimited by the two thin vertical lines in Fig. 1. In this 
mass range recoupling occurs for T <C m a as we can see in Fig. 1 where the diagonal dashed 
line denotes T = 3m a . Notice that here recoupling means axion decay since the putative ther- 
mal axion population is exponentially Boltzmann suppressed. Entropy is generally produced 
in the recoupling process, which in this case can be quite sizeable. In the approximate m a 
range 1-5 keV, axions decouple during the e + e~ annihilation epoch, but once more recouple 
in the Boltzmann suppressed regime. For m a = 5-20 keV axions recouple under relativistic 
conditions. 

For axion-photon interaction strength 5^1 the overall behaviour is similar. Both 
the Primakoff and decay rates are proportional to 5 2 and are increasing functions of m a . 
Therefore, a value 5 < 1 can be compensated by a increase in m a and vice versa, so the 
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Figure 1. Axion decoupling and recoupling (<5 = 1, C e = 1/6). Thick solid line: Freeze-out of 
Primakoff process. Medium solid line: Coupling and freeze-out of the Compton process. Thin solid 
line: Recoupling of inverse decay. In the yellow shaded region, axions are in thermal equilibrium, 
where the lighter yellow region is only relevant if the Compton process is effective. The dashed 
diagonal line denotes m a = 'ST. The vertical lines delimit the e + e~ annihilation epoch. 



decoupling and recoupling curves move upwards when 5 < 1 and downwards for 5 > 1. The 
same holds for the coupling to electrons C e which we have taken to be 1/6 in Fig. 1. 

3.1 Primakoff decoupling and axion abundance 

Next we must find the axion population existing after Primakoff /Compton decoupling in 
those cases when axions actually leave thermal equilibrium, i.e. when the cosmic evolution 
in Fig. 1 cuts through white space. We must solve the Boltzmann collision equation for the 
evolution of the axion number density n a during the e + e~ annihilation epoch. We cast the 
Boltzmann equation in the form 



dY a _r e + r ry ^ a dlogs 



dlogT H dlogT 3 



(Ya-Y^) , (3.5) 



where Y a = n a /s and s is the entropy density. We have numerically solved for the evolution 
of Y a and show the final axion abundance after e + e~ annihilation in Fig. 2. Hadronic axions 
with m a > 20 keV indeed keep thermal abundance. For smaller m a , the final abundance is 
smaller because axions decouple before the e + e~ entropy is fully released: they end up colder 
than photons (which define the equilibrium temperature). If axions couple to electrons the 
situation is very much the same except that the Compton process being more effective the 
abundances are generally larger than in the purely hadronic case. 

Axions modify the photon abundance in crucial ways whereas neutrinos are unaffected. 
Therefore, it proves convenient to normalize abundances in terms of the neutrino temperature 
which redshifts as T v oc i?" 1 after the neutrino decoupling. So we parameterize the axion 
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Figure 2. Axion number density n a after e + e~ annihilation from numerically solving the Boltzmann 
equation until T = m e /10. The equilibrium density n® q is defined in terms of the photon temperature. 
Solid line: only Primakoff process (hadronic axions). Dashed line: Primakoff and Compton. 

and photon abundances in terms of their temperatures as 

A =(!)'' and B=(g)\ ,3.6) 

As initial condition we use To = T v $ = T a o = 2 MeV when neutrinos have decoupled and 
e + e~ annihilation has not yet begun. The final epoch is defined by T\ = m e /10. Entropy 
conservation during e + e~ annihilation implies 

7 - + 2 + A = 2B X + A x , (3.7) 

where 7/2 and 2 are the e + e~ and 7 entropy degrees of freedom. From Fig. 1 we see that 
all of this happens far below the dashed line, i.e. axions are relativistic. We assume kinetic 
equilibrium even if axions decouple during the e + e~ annihilation. 

If axions decouple before e + e~ annihilation like neutrinos we have A\ = Aq. They do 
not receive any e + e~ entropy and photons are heated by the standard amount, i.e. B\ = 

On the other hand, if axions decouple during or after e + e~ annihilation they are some- 
what heated. We always have Aq = 1 and thus 

13 , A (T a y 

Bl = ^T2) h= B = \T) =^C=^- (3 - 8) 

It is h that was plotted in Fig. 2. If axions are fully coupled during e + e~ annihilation, then 
h\ = 1 and B\ = In general we will have -j- < B\ < 

3.2 Entropy production by axion decay 

Eventually axions disappear and transfer their entropy to the photon bath, increasing the 
photon density both relative to neutrinos and baryons. If axions never leave thermal equi- 
librium (for m a > 20 keV), we easily find that the final photon abundance is times the 
standard value. In general the entropy transfer to photons depends on two parameters, the 
initial axion abundance, parameterized by h, and the effectiveness of inverse decay when 
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axions become non-relativistic. To determine the final photon abundance we have solved 
numerically for the axion phase-space evolution (Appendix A). Here we only highlight the 
extreme cases of adiabatic decay (happening when inverse decay is effective at the time axions 
become non-relativistic) and out-of-equilibrium decay (when the inverse decay is ineffective) . 

If axions recouple relativistically, their temperature must catch up with photons and 
the latter are cooled. In this process, radiation is simply shuffled from one form to another 
and co-moving energy remains conserved. After axion recoupling, the photon abundance 
is reduced to B* = B\ [(2 + h^ 3 )/(2 + 1)] 3 ^ 4 . Comoving entropy increases by the factor 
[(2 + l)/(2 + hi)]B*/Bi, but this is never a big effect. For hi = 4/11 the entropy increases 
by some 2.6%. Later when axions get Boltzmann suppressed, this adiabatic process transfers 
their entropy to the photon bath, heating it according to -B2/-B* = (2+l)/2 = 3/2. Therefore, 
the final photon heating by axion recoupling and adiabatic decay is given by 



This ratio is the number density of photons relevant at CMB decoupling relative to the 
density they would have without axion effects, where in both cases the density is measured 
relative to neutrinos. Without axions B2/B1 = 1. 

If axions recouple non-relativistically, their out-of-equilibrium decay can produce a large 
amount of entropy. This point is illustrated by an analytic approximation to the entropy 
generation [35] 



where g*s is the effective number of entropy degrees of freedom and (gj s ) denotes an average 
over the decay time. In the last equality we have taken g*s = <7*s(3i) = 3.9. 

For our actual BBN and neutrino dilution arguments in the following sections we have 
numerically calculated the photon heating relative to neutrinos as a function of m a by solving 
the Boltzmann equation (Appendix A). In Fig. 3 we show the resulting BijB\ as a function 
of m a in the case 5 = 1. 

4 Big-bang nucleosynthesis 

The outcome of standard BBN depends basically on the parameter 77 = ns/n^, the baryon- 
to-photon ratio. The overall agreement of the predicted light-element abundances with ob- 
servations implies 5.1 < 7710 < 6.5 at 95% C.L. [29], where 7/10 = 10 10 r/. An independent 
determination t]\q = 6.23 ± 0.17 derives from the CMB temperature fluctuations [19]. The 
concordance between the BBN and CMB implied values will be disturbed in the decay- 
ing axion cosmology because the baryon abundance is diluted after BBN. In our scenario 
^bbn^cmb _ anc j therefore it follows the same curve depicted in Fig. (3). An in- 

creased baryon abundance relative to photons during BBN implies than nuclear reactions 
are more efficient, the deuterium bottleneck opens earlier, decreasing the residual deuterium 
abundance and increasing the yield of heavier elements. The increased expansion rate caused 
by the presence of axions works in the opposite direction but plays a sub-leading role. 

In order to quantify our arguments we have modified the publicly available BBN code 
PArthENoPE [36] to include the effects of axions, taking into account their impact on the 




(3.9) 




,l/3 )3 /4 




(3.10) 
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Figure 3. Photon density increase in our modified cosmology as expressed by B2/B1 for 6 — 1. Solid 
and dashed lines stand for hadronic and non-hadronic (C e = 1/6) axions respectively. The thin lines 
show the value if we assume that no entropy is generated in axion decay. 



Friedmann equation and the modified densities of different radiation species. Our results 
are shown in Fig. 4. For m a > 20 keV we treat axions as being in thermal equilibrium 
throughout BBN. For m a < 10 keV we use their abundance from our numerical freeze-out 
calculation, assuming that they are decoupled during BBN. We do not treat the intermediate 
case, leaving a gap in the predicted deuterium yield as a function of m a that is shown as a 
red band in Fig. 4. To calculate this curve we have adjusted the baryon abundance such that 
in the very end it matches the CMB-implied value. Its la range is represented by the width 
of the red band. 

We compare the predicted D yield with the measured value D/H| p = (2.82 ± 0.21) x 
10~ 5 [29], estimated from 7 high-redshift, low-metallicity clouds absorbing the light of back- 
ground quasars [37]. The dark (light) grey band in Fig. 4 is the experimental 95% (99%) C.L. 
region. We have included the independent determination of D/H| p from the 7 systems in the 
right panel of Fig. 4 as function of the H column density. Note that there is a significant 
scatter of the results whose origin still remains unclear and that has been taken into account 
by artificially increasing the error (the width of the grey band) [37]. Such a large scatter in 
the measurements of the D abundance may be a signal of some not yet understood processing 
of D inside these high-redshift clouds. Whether this is the case or not, there are no known 
astration processes that increase the D concentration, so the primordial should be larger 
than this estimation. As we see in Fig. 4, the presence of axions reduces ?/d, so our bounds 
appear to be conservative regarding this possible systematic. 

The deuterium abundance is reduced below its 2a observation if m a is below 300 keV. 
Therefore, BBN constrains axions to have masses 

m a > 300 keV . (4.1) 

Axions with mass above this limit have almost completely disappeared from the thermal bath 
before they can affect BBN and the predictions approach standard BBN. Since this bound 
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Figure 4. Deuterium yield yr> as a function of m a for S = 1 (Hadronic axions). The width of the 
red band represents the la uncertainty of the CMB determination of r\. To the right of the break, 
axions are treated as being in LTE, to the left they are assumed to be decoupled during BBN. The 
grey bands represent the 95% and 99% range for the observed D abundance. It is derived from 7 high 
rcdshift Ly-a clouds with individual results shown in the right panel as a function of the hydrogen 
column density [37]. 

corresponds to axions that are always in thermal equilibrium, even only via the Primakoff 
and inverse decay processes, it also applies to non-hadronic axions which would interact 
more strongly. Note finally that the helium yield is only mildly affected and does not provide 
interesting bounds. 

5 Neutrino dilution 

5.1 Radiation density at CMB decoupling 

We have seen that thermally excited axions that eventually disappear increase the abundance 
of photons at CMB decoupling relative to baryons and neutrinos. Since all cosmic parameters 
are defined relative to the observed CMB properties, this means that effectively the neutrino 
abundance is reduced. The impact of this effect on the CMB determination of the cosmic 
baryon abundance is minimal and thus could be neglected in the previous section. However, 
precision observables also measure the cosmic radiation density at decoupling, and here the 
"missing neutrinos" should make a difference. 

To extract cosmological information on the radiation density at CMB decoupling we 
analyse the usual 8-parameter standard ACDM model described in Ref. [8], extended in two 
ways. We allow the effective number of neutrino degrees of freedom to vary, assuming a 
flat prior on the interval < iY c g < 3.0. We recall that the radiation energy density is 
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Figure 5. 2D marginal 68%, 95% and 99% contours in the m v -N e g plane, where m„ is the individual 
neutrino mass (not the often- used sum over masses). 



traditionally expressed as 



Prad 



7 / 4 \ 4 ' 3 

1 + s{u) n <* 



Pi 



(5.1) 



where the photon radiation density is /3 7 = (7r 2 /15)T 4 and N e g the effective number of 
thermally excited neutrino degrees of freedom. The standard value is N c g = 3.046 instead 
of 3 because of residual neutrino heating by e + e~ annihilation [38] but, given the current 
experimental uncertainty, we neglect this tiny correction. In addition, we allow for neutrino 
masses, assuming a common value m u for all flavours and a flat prior on < s2^/s2 m < 1. The 
other parameters and their priors are identical with those provided in Ref. [8]. Moreover, we 
use the same set of cosmological data, i.e., the WMAP 7-year CMB measurements, the 7th 
data release of the Sloan Digital Sky Survey, and the Hubble constant from Hubble Space 
Telescope observations. 

Marginalizing over all parameters but m u and N e s we find the 2D credible regions shown 
in Fig. 5. Marginalizing in addition over m u we find the limits 



'2.70 at 68% C.L. 
N cS > \ 2.39 at 95% C.L. 

2.11 at 99% C.L. 



(5.2) 



These limits are very restrictive because on present evidence cosmology actually prefers extra 
radiation beyond N e g > 3 [17-21]. Taking this possibility seriously, we would need to worry 
about two novel ingredients: axions and radiation in some new form. 

The PLANCK satellite, currently taking CMB data, is expected to boost precision 
determinations of cosmological parameters. It will measure the cosmic radiation content at 
CMB decoupling with a precision of about AN c q = ±0.26 or better and thus will clearly 
decide if there is extra radiation in the universe [39] . If it finds convincing evidence for extra 
radiation, much in cosmology will have to be reconsidered besides our axion limits. 
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Figure 6. Radiation density during the CMB epoch (S = 1). Solid and dashed lines are for hadronic 
and non-hadronic axions (C e = 1/6). The 95% and 99% C.L. lower limits from Eq. (5.2) are shown 
as horizontal lines. 



5.2 Axion bounds 

The cosmic energy density in neutrinos is modified by the factor (T I f x /T^ td ) 4 between the 
axion and standard cosmology. The results of Sec. 3 imply that this ratio can be expressed 
in terms of the quantity B2, the modified (T/T u ) 3 value after axions have disappeared 

/ 11 \ 4/3 

AW=3(— ) . (5.3) 

We show the variation with m a for 6 = 1 in Fig. 6. At sufficiently high m a , the inverse decay 
process keeps thermal equilibrium during the decay. Even if axions decouple during e + e~ 
annihilation they will re-thermalize with photons before decaying. In this case, we have an 
analytical expression 



11 2(hi + 2) 2\ 



4/3 



For m a > 20 keV we have hi ~ 1 and N e g reaches asymptotically the minimum neutrino 
dilution 3 (11/13) 4 / 3 = 2.401. At much lar ger masses, m a ~ MeV, axions would disappear in 
LTE before neutrino decoupling, leaving no trace in cosmology. Therefore, at m a ~ MeV the 
value of N e R shown in Fig. 6 reaches a plateau at the the standard value N c r = 3. Comparing 
with the cosmological limits of Eq. (5.2) we see that even the minimum neutrino dilution is 
only barely allowed at 95% C.L. and therefore disfavoured, but not credibly excluded. 

For m a < 20 keV, axions decay increasingly out of equilibrium, creating entropy and 
reducing the final N e g further. Comparing the calculated N c g in Fig. 6 with the observational 
limits implies 

m a > 3 keV at 99% C.L. (5.5) 

The BBN limits reach to significantly larger masses, but the neutrino dilution limits are still 
nicely complementary. 
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6 CMB distortions 



The primordial plasma is optically thick before recombination, yet the photons produced in 
axion decay can produce spectral distortions in the CMB, measured by the FIRAS experiment 
to follow a black-body spectrum with 0(1O -4 ) precision [40]. 

After e + e~ annihilation, photon thermal equilibrium is maintained by Compton scat- 
tering, double Compton scattering and bremsstrahlung. Compton scattering is the fastest 
of these processes, but it can only provide kinetic equilibrium while conserving the photon 
number density. Double Compton (DC) scattering and bremsstrahlung (BR) are slower but 
change photon number and thus achieve thermal equilibrium. After a redshift of approxi- 
mately z ~ 3 x 10 6 , DC and BR no longer occur in equilibrium and any injection of photons 
leads to a non-zero pseudo degeneracy parameter, \x. The FIRAS data restrict its value to 

H<9x 10~ 5 . (6.1) 

Therefore, even a minute fraction of axion decays occurring after this redshift leads to a 
measurable distortion. 

Assuming that fi is always small, we can cast its evolution as [41], 

dp = djla _ ( J_ + J_\ ^ 62 ^ 
dt dt V*dc *br7 ' 

where £dc an d £br are the fi relaxation time scales due to the indicated processes (Ap- 
pendix A) and dpa/dt is the rate-of-change of \i generated by Compton scattering from 
photon injection or disappearance resulting from axion decay or inverse decay. It is related 
to the change in photon energy dp-y and photon number <in 7 by [41] 

dp a = J- ( - 4^1 V (6.3) 
P 2.14 V P 7 ™ 7 / 

assuming that /i is small. To obtain our bounds we have implemented Eq. (6.2) in a numerical 
code describing the evolution of axion decay (Appendix A). 

Our results are shown in Fig. 7, where we plot the final fx value as a function of m a and 
5 for hadronic axions. In the left panel we fixed 5 = 1. In this case we find 

m a > 8.7 keV at 95% C.L. (6.4) 

This bound is robust because \x is a steep function of m a . In particular it does not depend 
strongly on h\ and therefore is not expected to change significantly if we consider axions that 
couple to electrons, for which h\ will be slightly larger. 

The CMB distortion effect depends sensitively on the axion-photon interaction strength 
for 6 < 1 (right panel of Fig. 7). Generally the spectral distortions get larger for smaller 5 
at a given m a . For small 5, the decay happens later, when the photon distribution is less 
protected against distortions. 

For large 5 the final fi changes sign from negative to positive with increasing m a while 
for 5 < 0.1 n is always positive since axions decay non-relativistically injecting more energy 
than photon number. Of course, because of the sign change in [i there exist some fine-tuned 
cases where the final fx can be accidentally zero. 

Finally, we note that for photon injection occurring below z ~ 10 5 Compton scattering 
can no longer establish kinetic equilibrium, and the spectral distortion no longer takes the 
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Figure 7. Left: Photon fi parameter after axion decay (5 = 1). The observational upper bound 
is indicated. Right: Contours of log 10 \fj,\ in the (5-ra a -plane where black/red corresponds to posi- 
tive/negative values of fj, (The thick black line corresponds to the boundary fi = 0). The shaded area 
is excluded. 

shape of a pseudo-degeneracy parameter. Rather, for 10 3 < z < 10 5 it can be represented by 
a Compton-y parameter, and for z < 10 3 the decay would show up as a bump in the infrared 
background. Most important for our purpose here is that the bound on fj, is weaker than on 
other spectral distortions and therefore our bound is robust even for m a well below 1 keV. 

7 Conclusions 

Many cosmological particle-physics constraints were traditionally based on the idea that the 
successful predictions for the light-element abundances by BBN exclude additional radiation 
in the universe. However, over the past few years it has become apparent that both BBN 
and precision cosmology favour extra radiation on the level of around one effective neutrino 
family. While this preference is not yet strongly significant, it prevents one from deriving 
restrictive limits on possible additional forms of radiation. 

We have studied a scenario where BBN still provides useful constraints on low-mass 
particles. Axions in the sub-MeV mass range would produce photons late enough to dilute 
baryons after BBN relative to photons. Therefore, the baryon abundance established at CMB 
decoupling by cosmological precision data is smaller than at BBN, leading to an unacceptably 
lowered deuterium abundance. We have shown that this argument requires the mass of axions 
to exceed around 300 keV. Additional arguments derive from the cosmic radiation content 
at CMB decoupling (neutrinos would be diluted relative to photons) and spectral CMB 
distortions. These bounds depend sensibly on the two-photon coupling, which determines 
the axion decay, rather than on the couplings to leptons and hadrons. 

In the case of hadronic axions, together with the hot-dark matter limit, we thus find 
that cosmology alone excludes the m a range 0.7 eV-300 keV. Our constraints are complemen- 
tary to well-known stellar evolution limits and laboratory bounds, but based on completely 
different reasoning. 
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Figure 8. Primakoff (a) and Compton (b) processes on electrons and/or positrons and inverse two- 
photon decay (c). In the Compton case, there is also an u-channel diagram contributing. 
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A Axion out-of-equilibrium decay 
A.l Boltzmann equation 

We here compute the heating of the photon bath due to axion decay, with concomitant 
neutrino and baryon dilution and induced spectral fi distortion of the CMB. The relevant 
axion interaction processes are Primakoff and Compton scattering on electrons and positrons 
and two-photon decay and its inverse (Fig. 8). The Boltzmann equation for the axion phase 
space density f a {k a ,t) is 



dfa _ , R df a 
dt a dk a 



-£ - k a H^ = [C e (k a ,f a ,T) + C^K, f a ,T)} (f? - f a ) , (A.l) 



where H is the Hubble parameter. Here and in the following we denote boson four vectors 
as K = (cj, k) and k = |k|, whereas for fermions we use P = (E,p) and p = |p|. Specifically 
we use oj a = (k% + m 2 ) 1 / 2 . The collision integrals are 

C e = — [ , tg± , , ^7 (27r )^4 ( p + R p R ) \ M \2 

2oj a J (2vr) 3 2£i (2nf2E 2 (2vr) 3 2u; 7 V ; V 7 

x WBil/.Wtl - / 2 (£ 2 )]// a c >a) , (A.2) 

1 r (i 3 ki d 3 k 2 
2^7 (27r) 3 2wi (2vr) 3 2w 2 

x[l + /iW + / 2 N], (A.3) 



C 7 = 7T- / fo^sL (^Vl + K *~ K a )\M, 



where the particles are numbered as shown in Fig. 8. We assume thermal equilibrium between 
electrons and photons (temperature T), neglecting possible chemical potentials (see below). 
Relevant thermal effects are an effective photon mass, preventing (inverse) decay at high T 
and screening effects for the Primakoff process. The spin- and polarization-summed squared 
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matrix elements are [42] 



\M e \ 2 

\M P \ 2 



\M C \ 2 
2Re{M* P M c } 



\M,\ 2 



\M P \ 2 + \M C \ 2 + 2Re{M* P Mc} 



4:irag 2 



a i 



(t - m 2 ) ' 



-2m 2 mt 



2t 2 (s 



in 



ml) 



i + 2{s-m 2 e ) 2 -2m 2 a {s + m 2 e 



'}• 



2 ft 2 -m 2 a (3t-Am 4 e 



-16Trag ae g a -y 



(s — m 2 )(m 2 — u) 
{m 2 - t) 2 



m„ 



s + mi 



+ 



u + mi 



(s — Trig) 2 (u — m 2 ) 2 



(s — m 2 )(m 2 



1 



(A.4) 

(A.5) 
},(A.6) 
(A.7) 
(A.8) 



where we use the Mandelstam variables s = (Pi + Ky) 2 , t = (K^ — K a ) 2 and u = (K a — Pi) 2 . 
Also we have used g ae = (C e m e )/ f a . We keep the photon mass only in the propagator to 
account for screening effects and on the external legs in (inverse) decay to retain a threshold. 
For inverse decays we find analytically 



mi 



4m£ 77 



"7 i "-a 



mi, 



2T 1 
1 + — log - 

Pa 1 



-(w Q +p a )/2T 



-{u a - Pa )/2T 



a— >77 • 



(A.9) 



The Primakoff and Compton terms are much more complicated and have to be handled 
numerically. In this case it is much simpler to compute axion decoupling directly based on 
the Boltzmann equation for the axion number density, cf. Eq. (3.5). Integrating Eq. (A.l) 
over the axion phase space and defining Y a = n a j s (total entropy s) we find 



dYa 

' dt 



d 3 K 
(2vr) 3 



(c e + c 7 )(/ a -/ a eq ) = (r e + r 7 )K 



n. 



cq\ 



(A.IO) 



(A.11) 



where we have used comoving entropy conservation, d(sE?)/dt = 0, and assumed that f a is 
proportional to fa q , a reasonable approximation for computing the decoupling of a relativistic 
species. The second equality defines the rates we have used in Sec. 3. Entropy conservation 
can also be used to obtain the temperature-time relation that allows us to use the former as 
the independent variable 

d_ / d\ogs Y^ 1 d 

dt~ \dlogT 3 ) dlogT' 

which, together with Eq. (A.IO), leads to Eq. (3.5). 

Our bounds do not depend on a precise determination of hi, the axion abundance after 
the e 1 * 1 annihilation epoch, so we have solved Eq. (3.5) under two simplifying approximations. 
First, we note that the axion contribution to s is relatively small (the contribution from 
photons and neutrinos are at least a factor 4.4 bigger) so we can neglect it in s and in its 
derivative. Second, we have assumed Maxwell-Boltzmann distributions for all particles to 
find the average Primakoff/Compton rate, whereas decay and inverse decay are tractable 
without this simplification 2 . It is then easy to obtain numerically the evolution of Y a and 
the value of hi at T\ = m e /10. 



2 We have checked that our results do not differ significantly with the full results derived in [43, 44] for the 
relativistic limit T 3> m e , m a . 
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A. 2 Axion decay 

At T ~ m e /10 essentially all the e + e~ entropy has been transferred to photons and ax- 
ions. The e + e~ number density is exponentially suppressed, so we can ignore Primakoff and 
Compton. To compute numerically the photon heating when axions finally decay we evolve 
the axion phase-space distribution. The initial condition is set at T = m e /10 by computing 
the axion temperature that gives n a = hin^ when assuming f a to be a thermal distribution, 
albeit at a different temperature than photons. Note that the Boltzmann equation (A.l) 
can be much simplified by considering comoving momentum k a = k a (R/Ro), where Rq is a 
reference scale factor. Defining g a (k a ,t) = f(k a ,t), the l.h.s. of Eq. (A.l) simply is dg a /dt 
and the Boltzmann equation describing the decay becomes 

^ = -(fo-flS q )G y> (A.12) 

with C 7 given by Eq. (A. 9). 

We next compute the energy transfer to photons. The Boltzmann equation for the 
photon phase space distribution / 7 (/c 7 ,t) can be written as 

X [1 + /,(«,) + (/a" - /«) 

+ interaction terms with relic baryons and electrons. (A. 13) 

The assumption that photons are always in kinetic equilibrium during axion decay is equiva- 
lent to saying that the two sources on the r.h.s. will readjust themselves such that /~ is always 
of Bose-Einstein form. The evolution of the small photon chemical potential is considered 
below and neglected in the axion decay calculation. 

Multiplying by the photon energy cj 7 and integrating over phase space we obtain the 
equation for the evolution of the photon energy density 

^ - AHp, = J |^C 7 (g a - fl?) co, . (A.14) 

It states that the energy gain or loss of photons is due to axion decays or inverse decays. 
We have neglected the terms involving baryons and electrons. The evolution of the photon 
number density can be obtained analogously, but one has to take into account that double 
Compton and bremsstrahlung can yet modify the photon number (see below). 

To proceed further we define suitable time coordinates. As stated in the main text, 
neutrino temperature T v provides a natural reference to write the photon temperature in 
terms of T v through the parameter B that was defined in Eq. (3.6). Introducing the variables 
x = m a /T u and y = P/T u we have 

4- = xh4- , (A.15) 

dt dx ' v ' 

with H defined in Eq. (3.3). The evolution of g a in Eq. (A.12) is in final form 
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Furthermore, the comoving photon energy density p^R 4 = 7r 2 (Ti?) 4 /15 can be expressed in 
terms of B 4 . After some substitutions we find 



dx ~ ^] 2tt 2 H x {9a 9ah {AAV 

So far we have ignored the effects of the photon chemical potential. As explained in 
Sec. 6 the evolution of a small [i is given by Eq. (6.2). The characteristic time scales for 
double Compton scattering and bremsstrahlung to erase \i are given in Ref. [41]. Translated 
to our notation they are 

(rp \ 9/2 p3 -. 

. (A. 18) 

/ rp \ 13/4 D 9 / 2 , 

t BR = 3.4 x io 25 s (i - V2)" 1 (n b hr V2 {^f L ) ^ti^ti . 



with Y p the helium mass fraction and Q b h 2 the baryon density normalized to the critical 
density today (here h stands for the usual dimensionless Hubble constant). We use the values 
Y p = 0.23, r 7ito day = 2.725 K and tt b h 2 = 0.0223. Using Eq. (A. 14) and its counterpart for 
photon number, only retaining the axion decay induced dn 7 , the injection can be written as 
a time-dependent source for the pseudo degeneracy parameter 

dfia 1 f dyy 2 C 7 y^TF ( 45 Sir 2 

dx 2.14 J 2tt 2 [9a 9a) H x \^B* ^)B*y/W+tf 

We have numerically solved the evolution equations for the axion phase space distribu- 
tion Eq. (A. 16), the photon temperature Eq. (A. 17), and the pseudo degeneracy parameter 
Eq. (6.2), together with the initial conditions hq = and for B\ given by Eq. (3.8) as a 
function of h\. The latter was computed before, cf. Fig. (2). Note that the expressions for 
£dc and £br.) cf. Eq. (A. 18), require the value of -B to day> so in order to compute the chemical 
potential we first compute the evolution of B and determine -Btoday and then run again the 
evolution to compute the final value of p. 
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